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In this work we analyse the Parisi's oo-replica symmetry breaking solution of the Sherrington - 
Kirkpatrick model without external field using high order perturbative expansions. The predictions 
are compared with those obtained from the numerical solution of the oo-replica symmetry breaking 
equations which are solved using a new pseudo-spectral code which allows for very accurate results. 
With this methods we are able to get more insight into the analytical properties of the solutions. 
We are also able to determine numerically the end-point £ max of the plateau of q(x) and find that 
lim T -*o x max (T) > 0.5. 
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I. INTRODUCTION 



Since its proposal in the 80's the behaviour of the Parisi oo-replica symmetry breaking (oo-RSB) solution of the 
Sherrington-Kirkpatrick model has been extensively investigated both qualitatively and quantitatively jj], § . Despite 
this enormous amount of work, which has revealed many of the properties of the solutions, a complete control of the 
solution is still missing. One of the reasons can be traced back to the fact that till now only low order expansions were 
used, moreover applied often to reduced forms oo-replica symmetry breaking equations valid only near the critical 
temperature. From the numerical point of view there are only few works which confirm the general properties of the 
solution but do not allow for high accuracy. On the other hand oo-replica symmetry breaking solutions of the type 
encountered in the SK model have been found in other models of interest in different fields, e.g., in computer science 
with solvability problems || or in the study of the structural glass transition [0, pi @|. 

Motivated by these problems we believe that it would be quite useful to have some reliable and efficient tool to find 
good approximations of the full solution also far from the critical points. In this work we reconsider two approaches. 
The first one is based on expansions for temperatures near the critical temperature T c . As said above previous works 
considered only low order expansions |?], |[ Here, by using algebraic manipulators, we push the expansion to rather 
high orders and resumming it via Pade resummation technique we are able to a have good estimate of the solution 
for a wide range of temperature below T c . 

The second approach is numerical. Previous numerical studies of the oo-replica symmetry breaking solution used 
a naive integration scheme based on the direct discretization of the Parisi's equation [ll| [l3|. The main 
disadvantages of this approach are the large amount of memory needed for a good resolution of the solution and the 
numerical problems arising when q{x) is small. To overcome these problems we developed a new numerical scheme 
based on a pseudo-spectral algorithm which allows for rather accurate results for all temperatures with a reasonable 
amount of memory. Moreover the use of pseudo-spectral methods makes the whole code rather fast. 

We stress that while the methods we are going to discuss are applied here to the Sherrington - Kirkpatrick model, 
they have a wider range of application. In principle can be applied to any model with oo-replica symmetry breaking 
type solution j|. 

We find that for the Sherrington - Kirkpatrick model the Parisi solution q{x) is not an odd function as one may 
expect from its physical meaning. At any T < T c , the Taylor expansion of q{x) in powers of x contains both odd 
as well as even powers of x. The only term which is missing is x 2 . The presence of the fourth oder derivative was 
first noted by Temesvari Often, instead of q(x), it is more useful to consider the overlap probability distribution 
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function P(q), which gives the probability of finding two states with mutual overlap q according to the Gibbs measure. 
The two quantities are related by jl5L |l6) : 

where x{q) is the inverse function of q(x). In the absence of external magnetic fields the function P{q) must be an even 
function of q. The computed function q(x) is however defined only for positive values, therefore it determines only 
the right branch of the function P(q). If we define P(q) — dx/dq for q > then full P(q) is given by the symmetrized 
expression 

p(q) = \p{-<i) + \p{<i) (2) 

It is easy to see that the presence of non-zero even derivatives of q{x) at x — makes the function P{q) non analytical 
at q = 0: 

P(q) = c Q + c 2 q 2 + c 3 \q\ 3 + . . . (3) 

so that P(q) has discontinuous derivatives at q = 0. 

We shall discuss two different methods of computing the expansions. The first, discussed in Section y, performs 
expansion before imposing stationarity of the free energy functional. The two steps however can be inverted, i.e., the 



expansion can be done after stationarity is imposed, Section III. The two approaches are obviously equivalent and 



the advantage of using one or the other only depends on which quantity one is interested in. Since the expansions 



are likely to be asymptotic some resummation scheme, such as Pade discussed in Section IV, are needed. Finally 
in Section ^ we present a new integration procedure and compare the analytical results with those obtained from a 
direct numerical solution of the oo-replica symmetry breaking equations. 



II. EXPANSION OF THE FREE ENERGY FUNCTIONAL 

The Parisi's free energy / for a the SK model in an external field h at temperature T is [nj] : 



/= - (l-2o(l)+ / dxq 2 (x))+[ rxn I -- 

4 V q{ ' Jo q[ 'J J-oo V^Z(0) \ 2 g(0) 



djl 



(y - hf 



where 0(0, y) is the solution evaluated at x — of the the Parisi's equation 

q(x) 



with the boundary condition 



4"(x,y) + [3xcj ) '(x,yy 



(l,y) = /3- 1 log(2cosh/3y) 



(4) 



(5) 



(6) 



where we have used the standard notation and denote derivatives with respect to x by a dot and derivatives with 
respect to y by a prime. The order parameter q(x) at temperature T is obtained by the stationarity condition of (^) 
with respect to variations of q(x), while the value of at the stationarity point gives the free energy f(T). 

To expand the free energy functional (|4|) in powers of r = T c — T = 1 — T we observe that in the absence of external 
fields q(x) is different from q(l) only in a region [0,x max ] with x max = 0(t) |Q, so that an expansion in power of t 
must correspond to an expansion of the same order in x. Therefore to compute the free energy to order n we insert 
into eq. (0) the following expansions: 



q(l) = '£a i T i 



(7) 



and 



n—S n—3—i 



i=l j=0 



(8) 
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The coefficients of the expansion of the function <p(q, y) about q = q(l) and y = can be obtained by repeated 
differentiation with respect to q of the equation 



dq 



d 2 t (d<f> 

W + X{q) {dy- 



(9) 



Differentiating j times with respect to y this equation, mixed derivatives (q, y) can be eliminated in favor of 
derivatives with respect to y only. In the absence of an external magnetic field the last term in equation (Q) reduces 
to 4>(0, 0) greatly simplifying the calculation since at each step we can eliminate all terms containing odd derivatives 
of (j) with respect to y, as for example (d<p/dy) 2 in the previous equation, since all these vanish if evaluated at y = 
being 4>(q, y) and even function of y. 

Collecting all terms with the same power of r the free energy functional fl) is written as 



(10) 



This expression must be stationary with respect to variations of a's and 6's for any r. Imposing stationarity of each 
Ci we can find the value of the parameters a and b. For example to order t 6 we have: 



_I + ^ t + 3 t 2 + 38t j I , 



+ (-1-9T-30- 



351 9189 \ 5 27 6 



and 



Zmax = 2t - 4t 2 + 12t 3 - 69 t 4 



2493 - 20544 , 

T 5 T b . 

5 5 



(12) 



By using this procedure we have obtained the free energy up to order 30, q{x) to order 13 and q(l) to order 14 because 
despite the fact that the free energy is evaluated to order n, the variational relations allow to determine x(q) only to 
order [(n — 3)/2j and q(l) only to order [{n — l)/2]. 

From eq. (O) we clearly see that q(x) contains even powers of x, with the exclusion of x 2 . In the next Section 
we shall derive exact relations among the derivatives of q(x) at x = from which follow that q^ 2 \x = 0) = but 
q^{x = 0) i= 0. 



III. EXPANSION OF THE ORDER PARAMETER q{x) 



To evaluate the derivatives of the order parameter q(x) at x = we use a variational approach developed by 
Sommers and Dupont|ll[]. This method as also the advantage of leading to exact relations among derivatives of 
different order, so can be used to test the findings of the previous Section in a non-perturbative way. The starting 
point is the variational form of the Parisi's free energy /: 



/ 



|(l-2g(l)+ / ./..■,,-(,•)) + 







dy 



V2^(0) 



exp 



(y - hf 

2<z(0) 



my) 



dyP(l,y) [0(l,y)-Tlog(2cosh/3y)] 

— oo 
1 r+oo 

dx / dy P(x, y) 

J -oo 



<l>(x,y) 



q(x) 



cf>"(x,y)+Pxcb'(x,yf 



(13) 



Imposing stationarity with respect to variations of P{x, y), P(l, y), <^>(a;, y), <j>(0, y) and q{x), one obtains the variational 
equations: 



m(x,y) 



q{x) = J dy P(x, y) m 2 (x, y) 
q(x) 



m! (x, y) + 2 j3xm(x, y) m'(x, y) 



(14) 
(15) 
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P(x, y) = ^ [P"(x, y)-2(3x [m(x, y) P(x, y)}'] (16) 

with initial conditions (in the absence of a magnetic field) 

m{l,y) = tanh(y/T) (17) 
P(0,y) = S(y) (18) 

These equations are the starting point of both the expansion discussed in this Section and the numerical solution. 

A time scale t x can be associated to the order parameter q(x) such that for times of order t x states with an overlap 
equal to q(x) or greater can be reached by the system. In this picture the P{x, y) and m(x, y) become respectively the 
probability distribution of frozen local fields y and the local magnetization in a local field y at the time scale labeled 

The derivatives of q(x) can be obtained by successive x-derivation of eq. (|14|). The procedure is simplified by the 
use of the following identity jt| : 

— J dy P(x, y) f(x, y) = J dy P(x, y) Sl(x, y) f(x, y) (19) 

where 

n ^-T x + \{w + 2Pxm(x ^y) m 

The application of the operator £l(x, y) generates derivatives of the function m( Xl y) with respect to x and y. Mixed 
derivatives such as rrS 1 ^ (x,y) can be eliminated in favor of derivatives of m(x,y) respect only to y by deriving 
equation ( [l5| ) j times with respect to y. 

The first application of this procedure yields 

q{x)=q(x) J dyF(x,y)(m') 2 (21) 

which for q(x) =/= simply reads 

1 = J dyP(x 1 y)m'( X ,y) 2 . (22) 
The procedure can be iterated infinitely. For example, the next three applications leads respectively to 

= ~ f dyP{m'f + f dyP(m") 2 (23) 



and 



Pdy 



dyP{m'f = qj Pdy (V") 2 - ^m'(m") 2 + 6 [^f (rn') 4 ) (24) 
f(18xq + 6x 2 q)(m') 4 (18<?T + Ylx'qT - 120m'x 2 q 2 )m'(m") 2 



\ T 2 



-W X q 2 {m") 2 +qm'"T-20xq 2 m l {m l ") „, 24x 3 g 2 (m') 5 , 2 2 , 

-m — + q (m ) =0 (25) 



T T 3 

We are interested into the derivatives of q(x) at x = 0, so we take the limit x — > of the above equations. The limit 
can be done in trivial way and, since the function P(0,y) reduce to a (S-function [see eq. (Esl)], the equations are 
greatly simplified. Moreover since m(x, y) is an odd function of y for any x clearly ^> (0, 0) = for any even j . In 



this limit equations (|22|), ( |24[ ) and fl25[ ) reduce respectively to: 

1 = m'(0,0) (26) 

|m'(0,0) 3 = g(0)m"'(0,0) 2 (27) 

g(0)m"'(0,0) 2 = (28) 
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while equations ( [14 ) a nd (E3)become trivial identities. 
From equations (26|) and (27J) we have 



m'"(0,0) = -^^0 (29) 

therefore ( p8| ) implies that g(0) = as already found in Ref. ||. 

To obtain information on the fourth derivative of q(x) the above procedure must be iterated two more times. Since 
successive derivatives yields expressions with a rapidly growing number of terms we only report the limit x — > result: 

^^^(O.O)' - 384(0) 2 ^(0,0) 2 3 = Q 

T 1 T 



(4)(oK , (0;0) J6,(QK(0,0)3^ q (3i) 



where equation ( p6|) and the exact result g(0) = have been used. Note that equation ([31]) , with equation (|29|), gives 
a complete determination of the quartic derivative of q(x) at x = as a function of the temperature T and of the first 
derivative q(x = 0): 

= (32) 

This relation shows that the function does not have a well defined parity |l4j ]. 

Going to higher orders one can show that all the even derivatives can be expressed in terms of the odd ones; for 
instance we have 

(6) 34272V2g(0) 7 / 2 1056V2g(0) 3 /y 3 > (0) 

q W _ T 5/2 T 3/2 

and so on. 

In the limit T -> we have Tg(0) = 0.743 ± 0.002. Note that if we take q(0) ~ 1/T for T -> the previous 
equations implies that all the derivatives diverge with the temperature as q( n \Q) ~ l/T™, in agreement with the 
Parisi- Toulouse scaling q(x,T) — q(f3x) jl^, |T^]. Note that we have derived this scaling under strong hypothesis that 
it must be valid asymptotically for T — > and fix — > 0. 

This approach also provides an alternative method to compute the expansion of q(x) in powers of x and t: starting 
from q(x) evaluated at a given order in x and r we can compute m^'^O, 0) through ( |l5| ) and then g(a;) at the next 
order through the set of equations (|2^),(^), (^),(|3^) and so on. The set of equations can be solved iteratively. By 
this method we were able to compute the series expansion of q{x) up to order 20, improving the results of previous 
section. 

IV. RESUMMATION OF THE EXPANSIONS 

Unfortunately all the expansions derived in the previous Sections are likely to be asymptotic and to obtain sensible 
estimates of the various quantities of interest some form of resummation must be done. Here we shall consider the 
standard Pade approximants which for a series of degree N + M reads |19] : 



■\7V 

l + £f=i 



P£W - : .. (34) 



where the coefficients are chosen so that the first (N + M + 1) terms of the Taylor expansion of P^(x) match the the 
first (N + M + 1) terms of the of the original series. In the following we shall call this the Pade approximant (N, M). 

The first problem we faced is that despite the fact that the series have alternate signs, they are not Stijlties integral 
and therefore we cannot obtain in a systematic way a sequence of lower and upper bounds p"9] ] . This difficulty can be 
overcome by noticing that most of the quantities we are interested in, such as for example free energy or entropy or 
q(l), do have a null derivative at T = 0. Therefore an indication on the quality of the approximants can be obtained 
by analyzing the behaviour near T = 0. For example, the free energy as a function of T is reproduced quite well by 
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-0.763 



-0.7631 - 




FIG. 1: Free energy as function of r = 1 
(10,11) and (17,10). 



T for different Pade approximants. Top to bottom: (14,16), (13,12), (12,11) 



many Pade approximants, even at very low orders, however some of these have a positive derivative at T = while 
others negative, see Fig. [l]. By inspecting the figure we can safely assume that approximants with pos itive derivative 
give an upper bound, while those with negative derivative a lower bound, for the true free energy f20fl . 

As a general fact we obtain that the best Pade approximants at a given order in r are those with nearly the 
same degree of the numerator and the denominator. We stress, however, that as usual with asymptotic expansion an 
increase of the order in r does not necessarily correspond to an improvement of the precision. With this procedure 
we obtain for the free energy an estimate with at least 16 digits precision at T = 0.9 and 8 digits at T = 0.5, and 
for the ground state energy E = —.76321 ± .00003 in agreement with Parisi's estimate E = —.7633 ± .0001 [ JPtJ . A 
similar analysis can be used to determine the value of x max as a function of temperature, the result is shown in Fig. 
H. The value of the breaking point is finite in the limit T — > 



c (0) = .548 ± .005, 



(35) 



see inset Figure |2j, and slightly larger than the value 1/2 predicted by the Parisi- Toulouse scaling, in agreement with 
the approximate nature of this relation 



K 
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The analysis of the function q(x,T) is more complex, because not only the Taylor expansion of q(x) in powers 
of x around any < x < x max is likely to be asymptotic for any fixed temperature, but the expansion in r of the 
coefficients of the x-expansion are themselves non convergent. Therefore one should use a double Pade expansion, one 
for the coefficients and one for q(x). The procedure however is quite difficult because we do not have a systematic 
way of choosing the best approximant and, moreover, coefficients of higher order are known with less precision in r. 
A better approach is to construct the function q(x) directly point by point by computing q(mx max ) where m = i/n, 
(i = 0, 1, . . . n) for fixed n. For any m and T the quantity g(mi max ) is itself a power series in r which can be summed 
up using Pade approximants. With this procedure the function q{x) can be determined for different x-resolution just 
changing the value of n, e.g., n = 50, 100, 1000, and using the value of x max previously found, see Fig. ^. In Figure || 
the function q(x) is shown for various temperatures T . 

This method can be extended to any function of x or q, for example, we computed the overlap probability function 
P(mq max ) in a wide range of temperature T > 0.3, see Figure ^. We found that the best Pade approximant is given by 
the (12, 7). By using the relation (|IJ) we can have an independent estimation of q(x) with which to test the precision of 
our results. By using a norm dca(q, q') — maxo<a;<i \q(x) — q'{x)\ and expansions up to order 20 we find, for example, 
that doo(g,9') = O(10- 5 ) for T = 0.6 and d^q') = O(10- 4 ) for T = 0.4. 

The form of the function q(x) confirms the prediction of Ref. |l(J obtained from interpolation of the 11-RSB 
solution. In particular it confirms the approximate scaling q(x.T) ~ q(x/T) at low temperatures, see Figure^. Note 
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that the scaling fails when (3x ~ O(l), in agreement with the findings of previous section. 

Finally we mention that an alternative resummation technique based on the Borel transform give results consistent 
with those obtained with the Pade approximants. 




FIG. 4: P(q) as function of q for various temperatures. From left to right T = 0.95 to T = 0.30 in step of 0.05. The data are 
obtained with a (12, 7) approximant. Note that P(q) attains its minimum value for q > 0. This happens for any temperature 
T < 0.961938.... 



V. NUMERICAL INTEGRATION OF THE oo-RSB EQUATIONS 

To check the analytical results of the previous sections we have solved numerically the oo-RSB equations ([TI]) - jilil ) 
on a discrete set of points in the infinite strip [0 < x < 1; — oo < y < oo] and determined q(x), P(x,y) and m(x,y). 
The numerical method is based on the iterative procedure of Ref. fl2|| : from an initial guess of q(x) the fields m(x, y), 
P{x, y) and the associated q{x) are computed in order as: 

1. Compute m{x 1 y) integrating from x = Xo to x = eqs. ( |i~5| ) with initial condition (|l7|). 

2. Compute P(x,y) integrating from x — to x — x eqs. ( |l6| ) with initial condition (|i"§|). 



3. Compute q{x) using eq. (14) 



where xo < 1 (See later). The steps 1. — > 2. — > 3. are repeated until a reasonable convergence is reached, typically 
mean square error on q, P and m of the order O(10 -6 ). 

The core of the numerical scheme is the integration of the partial differential equations ( |T5| ) and ( [l6| ) along the x 
direction which, at difference with previous numerical studies |l^, , is done in the Fourier Space of the y variables 
where the equations take the form: 



d k 2 q{x) /3q(x) r 2 -, 

m{x, k) — — - — m{x, k) — \k FT [m \ [x, k) (36) 



dx 
and 

^P{x,k) = - k ^ P(x, k) - (3q(x) ifc TT [P m] {x,k) (37) 

For each wave- vector k these are ordinary differential equations which can be integrated using standard methods. To 
avoid the time consuming calculation of the convolutions in the non-linear term we use a pseudo-spectral |^l| code on 
a grid mesh of N x x N y points, which covers the x-interval [0, xq] and the y-interval [— y max , J/max]- The truncation of 
wave-number may introduce anisotropic effects for large k, therefore to ensure a better isotropy of numerical treatment 
we perform de-aliasing via a N y /2 truncation [ p2| . Finally the x integration has been performed using an third-order 
Adam-Bashfort scheme which has the advantage of reducing the number of Fast Fourier calls (23]] . Typical values used 
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q 




FIG. 5: q as a function of fix for different values of T. Top to bottom T = 0.30, T = 0.35, T = 0.40, T = 0.45 and T = 0.50. 

are N x = 100 5000, N y = 512 4096 and y max = 12 -j- 48. The difference between the values used for A^ and A^ 
follows from the observation that if the solution in the y-direction is smooth enough, then only few low wave-vectors 
are exited. The value of the parameter y max fixes the y-range where the solution is assumed different form zero, since 
in the numerical algorithm is implicitly assumed that 

P(x, y) = m(x, y) = \y\ > y max . (38) 

This explain the rather large value used. The number of iterations necessary to reach a mean square error on q, P 
and m of order O(10 -6 ) depends on the initial guess of q(x) but it is typically of few hundreds. 

In Figure ^| are shown the order parameter q(x) and the overlap probability distribution function P(q) at T = 
0.6 computed for increasing ^-resolution and xo — 1. As expected the agreement between the numerical and the 
perturbative solutions increases with the number of N x of x-grid points. However, the convergence is not uniform: it 
is rather fast far from x max and much slower for x ~ x max , see the inset of Figure ^ panel (a). This is not unexpected 
because for x — x max the derivative of the order parameter q^(x) has a cusp: 

lim g (1) (a;)>0, lim q {1) (x) = (39) 

X — *X max X — ^max 

making the convergence more difficult. We recall that in deriving equations ([l5]) and ([l6|) differentiability of q(x) was 
assumed. The use of lower order integration schemes, as second-order Adam-Bashfort or Euler schemes, does not give 
sensible improvements. 

Larger values of N x requires larger needs of computer memory therefore to increase the precision we adopted a 
different approach. Since q(x) — for x > x max equations (|l5|) and ([l6]) are trivial in this range and we can reduce 
the upper bound of the x-integration from x — 1 to x — xq — x max . This obviously requires the knowledge of x max 
for the given temperature. However if we assume no a priori knowledge of x max we must proceeds for successive 
approximations: we start from x$ = 1 an then reducing it until we 'hit' the value of a: max . This procedure is simplified 
by the fact that if xq < x max the shape of q(x) near xq changes dramatically with the concavity passing from negative 
values for xq > a; max to positive values for xq < x max . In Figure panel (a) are shown q{x) (panel a) and P(q) (panel 
b) at T = 0.6 for different values of xq, the improvement in rather evident. As additional check we have considered 
the equality 



1 — / dx q(x) = T 
Jo 



(40) 
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FIG. 6: q(x) as function of x (panel a) and P(q) as function of x (panel b) at T = 0.6 for (a: bottom to top, b: left to right) 
N x = 50, N x = 500, N x = 1000 and N x = 5000. In all cases x = 1, i/max = 48 and N y = 4096. The dash line is the result 
from the perturbative solution discussed in the previous sections. Inset panel a: enlargement of the region near x max . 



which is satisfied by our numerical solution for all studied temperatures with at least four digits. For example we for 
T = 0.8 we get 0.79999(4), while for T = 0.5 the value is 0.49999(3). 

Note that not only by fine tuning of xq we can have a good solution for q(x) at the given temperature, but we also 
have the value of xq. This is best seen by analyzing the concavity of q(x) near xq. In Figure ^ we show the second 
derivative of q(x) near xq for T = 0.4 and N x = 500, from which one may conclude that 0.505 < x max < 0.510. 

A careful analysis of the stability of this results as function of N x , see Figure ^, reveals, however, that the correct 
estimation is 0.510 < x max < 0.515, in rather good agreement with the analytical result x max = 0.5111 ± 0.0002. The 
same analysis for T = 0.6 leads to 0.438 < x max < 0.440. 

We are now in the position of checking the results of previous section about the derivative of the order parameter 
at x — 0, and in particular the conclusion 



lim q (3) {x) > 0. 



(41) 



In Figure |l^ we show the second and third derivative of q(x) obtained from numerical differentiation of q(x). The 
agreement with the per turbative result is sufficiently good, moreover from the right panel of Figure [l0| we clearly see 
that the prediction (|l]) is verified. 

We conclude this Section with a short discussion on the entropy which, using the stationarity of the free energy 
functional (IL3[), can be written as: 



/? 2 



dy P(l,y)[log2cosh/?2/-y tanh(/3y)] . 



(42) 



For other equivalent forms see, e.g., Ref. 
Figure 
find 
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. The entropy as function of temperature is shown in the left panel of 
The entropy must vanish quadratically with the temperature as T — > . From our numerical data we 



lim -V = a- 0.72 



(43) 



to be compared with 0.718 ± 0.004 of the analytic expansions. 
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FIG. 7: Panel a: q(x) for T = 0.6 near x ma x for different xo- Panel b: P(q) for T = 0.6 for different xo ■ In all cases N x = 500, 
2/max = 48 and N y = 4096. The dash line is the result from the perturbative solution discussed in the previous sections. 



In the limit T — ► the quantity 1 — q(l) must also vanish as T 2 | fll| . The behaviour of g(l) as function of T is 
shown in the right panel of Figure [ll]. Using this data we obtain 

li m 1 ~ y ) ~ 1.60 (44) 
T^O T 2 v 7 

in very good agreement the value 1.60 ± 0.01 obtained with the expansions of previous sections. 



VI. CONCLUSIONS 



In this paper we have studied the properties of the oo-replica symmetry breaking solution of the Sherrington - 
Kirkpatrick model without external fields. Using high order expansions in t = T c — T we are able to compute the 
order parameter q(x) and other relevant quantities for a large range of temperatures with high precision. In particular 
we found that q(x) is not an odd function of x, confirming the prediction of Ref. jl4j]. Direct consequence of this is 
that the overlap probability distribution function P{q) has discontinuous derivatives at q = 0. Another consequence 
of our findings is that the Parisi- Toulouse scaling becomes exact asymptotically for T — > and (3x — > 0, while for 
T — > is a fairly good approximation. This is also consistent with the T = limit of the breaking point which we 
found to be x max (0) = -548 ± .005. 

Having reached very high orders we can reasonably speculate on the analytical properties of the function ^(ar). In 
particular we believe that 

• All the expansions in power of r are asymptotic expansions. 

• At any temperature, the function q(x) is infinitely diffcrcntiable but not analytical for any x, in particular the 
Taylor expansion of the function q(x) around any < x < x max does not converge but is asymptotic. 

This singular behaviour is not connected neither with the replica limit nor with the Parisi Ansatz, actually it originates 
from the singularities in the complex plain of the initial condition of the Parisi equation: f(l,y) — In 2 cosh/3?/. This 
is clearly seen for the replica-symmetric solution 

/+oo ^ z ^ 
—=e^ tanh 2 : (/3 '^/qz) 
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FIG. 8: Second derivative of q(x) for T — 0.4 for different xo and N x 
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In this case it is easy to prove that the expansion of (1 — T 2 ) in powers of p = (3 2 q is asymptotic because it corresponds 
to substitute tanh 2 z in the integrand with its Taylor expansion which is not convergent on the whole real axes. Then 
one can prove that the expansion of q in powers of r = 1 — T is asymptotic recalling that standard manipulation 
(e.g. multiplication, division, inversion...) on an asymptotic expansion in power series do not change its character. 
A detailed treatment of the RSB solution is much more complex, but the origin of the asymptotic character is likely 
to be the same. Indeed an expansion in small r (and therefore in small q) corresponds to an expansion in small y of 
all the quantities like f(x,y) and m(x,y); the appearance of integrals of the form J Pfdywheve P(x,y) ~ exp(y 2 /a:) 
generates asymptotic expansions since the Taylor expansions of f{x, y) and m(x, y) in powers of y do not converge on 
the whole real axes. These arguments can be very useful in practice to guess the position of the singularities of the 
Borel transform if one want to sum the expansions through a conformal mapping [B4J . For instance in the expression 
of the free energy appear integrals of the following form: 



dz 



■.e 2t lncosh(z) 



(45) 



The singularities of the Borel transform of the previous integral are located on a cut running from — oo to — 7r 2 /8 and 
a possible guess is that this be also the singularity structure of the Borel transform of the free energy. This guess is 
supported by the analysis of the series expansions. 

The analytical results have been compared with numerical solutions of the oo-replica symmetry breaking equations. 
We have developed a new numerical approach based on a pseudo-spectral code which leads to a strong enhancement of 
the quality of the numerical results. We have also shown how, for example, to determine the value of x max numerically. 
In all cases the agreement between the numerical and the analytical results is rather good. 

We conclude by stressing that our results go beyond the interest on the Sherrington - Kirkpatrick model, since the 
method we used here are far more general and can be employed to a wider class of models with generalized oo-replica 
symmetry breaking equations such as those introduced in Ref. ||. In particular in this reference the numerical 
method was applied to the 3-SAT model, and the extension to other relevant models is under development. 



[1] M. Mezard, G. Parisi, M.Virasoro, Spin Glass Theory and Beyond, World Scientific, Singapore (1987). 



13 




0.52 



FIG. 9: Second derivative of q(x) for T = 0.4 for different N x and x = 0.510, j/ max = 48 and N v = 4096. Left panel: 2:0 = 0.510, 
increasing N x leads to a positive value of q^(xo) implying xq = 0.510 < :r max . Right panel: xq = 0.515, increasing N x leads 
to a more negative value of q ( - 2 \xo) implying xo = 0.515 > a; max . 

[2] K.H. Fischer and J. A. Hertz, Spin-Glasses (Cambridge University Press, 1991) 

[3] A. Crisanti, L. Leuzzi and G. Parisi, J. Phys. A: Math. Gen. 35, 481 (2002). 

[4] E. Gardner, Nucl. Phys. B257, 747 (1985); 

[5] T.R. Kirkpatrick and D. Thirumalai, Phys. Rev. B 36, 5388 (1987). 

[6] M. Sellitto, M. Nicodemi and J.J. Arenzon J. Phys. I France 7, 945 (1997). 

[7] G. Parisi, Phys. Rev. Lett. 43, 1754 (1979); J. Phys. A 13, L115 (1980). 

[8] I. Kondor, J. Phys. A 16, L127 (1983). 

[9] H.-J. Sommers, J. Physique Lett. 46, L-779 (1985). 

[10] J. Vannimenus, G. Toulouse and G. Parisi J. Physique 42, 565 (1981). 

[11] H. J. Sommers, W. Dupont, J. Phys. C 17 (1984) 5785-5793. 

[12] K. Nemoto, J. Phys. C 20, 1325 (1987). 

[13] P. Biscari, J. Phys. A 23, 3861 (1990) 

[14] T. Temesvari, J. Phys. A 22, L1025 (1989). 

[15] C. De Dominicis and A.P. Young, J. Phys. A 16, 2063 (1983); 

[16] G. Parisi, Phys. Rev. Lett. 50, 1946 (1983). 

[17] G. Parisi, J. Phys. A 13, L115 (1980). 

[18] G. Parisi and Toulouse, J. Physisque Lett. 41, L-361 (1980). 

[19] CM. Bender and S.A. Orszag, Advanced Mathematical Methods for Scientists and Engineers, McGraw-Hill (1978). 

[20] We note that it is possible to impose a null derivative at T = directly into the Pade approximants. This however does 

not produce a sensible improvement of the precisions at a given order, making at the same time more difficult to have a 

control on the convergence. 

[21] S. A. Orszag, Studies in applied mathematics (Cambridge University, Cambridge, 1971), Vol. 4, p 293. 

[22] G. S. Patterson and S. A. Orszag, Phys. Fluids 14, 2538 (1971) 

[23] See, e.g., J.H. Ferziger and M. Peric, Computational Methods for Fluid Dynamics, Springer- Verlag (1996). 

[24] G. Parisi, Statistical Field Theory, Addison Wesley (1988) 




FIG. 11: Left panel: entropy s as function of temperature T. Right panel: g(l) as function of temperature T. 



